%
%  This MATLAB/OCTAVE script reads Fourier components of a 
%  a field from SpectralPlasmaSolver (SPS) output, performs 
%  Fourier transform, and save the resulting representation 
%  of the field in physical space into a new file.
% 

% close all all open files and clear variables
clear all;close all

% quantities to read
q = {'bx','by','bz','ex','ey','ez'};

% time step to process
step = 1216

% where the SPS output is located
in_dir = 'data'

% where to save output 
out_dir = 'data-physical'

% spatial resolution of the simulation (number of Fourier modes)
Nx = 399;
Ny = 399;
Nz = 249;

% total number of modes
Nk = (Nx)*(Ny)*(Nz);

% loop over quantities
for iq = 1:length(q)

  % read input data: array of double-precision, complex numbers stored in
  % Fortran (column-major) ordering, in the format 
  % (real_part_1),(imaginary_part_1),(real_part_2)... 
  fname = sprintf('%s/%s_%i.bin',in_dir,q{iq},step)
  f = fopen(fname,'rb');
  tmp = fread(f,2*Nk,'double');
  fclose(f);
 
  % convert real array tmp into a complex, appropriately shaped array 
  field = tmp(1:2:end-1) + 1i*tmp(2:2:end);
  field = reshape(field,Nx,Ny,Nz);
 
  % take 3D Fourier transform. Note that Forier components are indexed 
  % as -(Nx-1)/2..(Nx-1)/2, so they must be shifted before applying ifftn
  real_field = Nk*real(ifftn(ifftshift(field)));
  
  % save the physical space data in single precision 
  fname = sprintf('%s/%s_%i.gda',out_dir,q{iq},step);
  f = fopen(fname,'wb');
  fwrite(f,real_field,'single');
  fclose(f);

end

